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Abstract We investigate how the divergence-free property of magnetic fields 
can be exploited to resolve the azimuthal ambiguity present in solar vector 
magnetogram data, by using line-of-sight and horizontal heliographic derivative 
information as approximated from discrete measurements. Using synthetic data 
we test several methods that each make different assumptions about how the 
divergence-free property can be used to resolve the ambiguity. We find that the 
most robust algorithm involves the minimisation of the absolute value of the 
divergence summed over the entire field of view. Away from disk centre this 
method requires the sign and magnitude of the line-of-sight derivatives of all 
three components of the magnetic field vector. 
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1. Introduction 

Measurements of the vector magnetic field in the solar atmosphere are cru- 
cial for understanding a wide variety of physical phenomena. Using the linear 
polarisation of magnetically sensitive spectral lines to infer the component of 
the field transverse to the line-of-sight results in an ambiguity of 180° in its 
direction (Harvey, 1969). This ambiguity must be resolved in order to completely 
determine the magnetic field vector. Several methods are currently in use to 
resolve the ambiguity (for overviews see Metcalf et al, 2006; Leka et ai, 2009). 
In general, each method involves an assumption or approximation that may not 
be appropriate for observations of solar magnetic fields. 
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A promising approach, which can potentially avoid unrealistic assumptions, is 
to exploit the divergence-free property of magnetic fields (e.g., Wu and Ai, 1990; 
Cuperman, Li, and Semcl, 1993; Li, Cuperman, and Semel, 1993; Boulmczaoud 
and Amari, 1999; Li, Amari, and Fan, 2007; Crouch and Barnes, 2008). One 
challenge in calculating the divergence is that it requires information about the 
variation of the field in the direction perpendicular to the solar surface. To this 
end it is possible to derive from observations the variation of the field along the 
line-of-sight direction (e.g., Ruiz Cobo and del Toro Iniesta, 1992; Collados et 
ai, 1994; Metcalf et ai, 1995; Del Toro Iniesta and Ruiz Cobo, 1996; Liu et ai, 
1996; Westendorp Plaza et ai, 1998, 2001; Socas-Navarro, Trujillo Bueno, and 
Ruiz Cobo, 2000; Eibe et ai, 2002; Leka and Metcalf, 2003; Socas-Navarro, 2005, 
2007). The line-of-sight direction is not perpendicular to the solar surface, except 
at disk centre; therefore some care is required to correctly treat the geometrical 
perspective. 

In Crouch and Barnes (2008, henceforth paper I) we showed that the ambi- 
guity can be resolved with line-of-sight and horizontal hcliographic derivative 
information by using the divergence-free property without additional assump- 
tions about the nature of the field. To resolve the ambiguity away from disk 
centre wc demonstrated that it is sufficient to determine four pieces of infor- 
mation: i) the sign of the line-of-sight derivative of the line-of-sight component 
of the magnetic field; ii) the sign and magnitude of the line-of-sight derivative 
of the magnitude of the transverse component of the field; Hi) the sign and 
magnitude of the line-of-sight derivative of the azimuthal angle; and iv) the sign 
and magnitude of the horizontal heliographic derivatives of the magnitude of 
the transverse component of the field and the azimuthal angle. Paper I was a 
theoretical examination in which we assumed that all necessary derivatives could 
be determined exactly at a given location and that the solar surface could be 
represented by the hcliographic plane. 

The aim of this paper is to investigate how the divergence-free property can 
be used to resolve the ambiguity for data where the magnetic field is measured 
at discrete spatial locations. Using synthetic data we test two methods that have 
previously appeared in the literature (e.g., Wu and Ai, 1990; Cuperman, Li, and 
Semel, 1993; Li, Cuperman, and Semel, 1993; Boulmezaoud and Amari, 1999; 
Li, Amari, and Fan, 2007); both of these are based on the divergence but they 
make different assumptions about how to use it to resolve the ambiguity. We 
show that both of these methods can produce significant errors. For this reason, 
we present a more robust algorithm, the global minimisation method. 

The outline of this paper is as follows. In Section 2 we describe the derivation 
of the divergence-free condition for any position on the solar disk in terms of 
observable quantities and discuss the consequences for ambiguity resolution. In 
Section 3 we review the synthetic data and metrics that will be used to test 
the performance of the ambiguity resolution algorithms. In Section 4 we test 
the algorithms and the validity of the assumptions made. In Section 5 we draw 
conclusions. 
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2. The Divergence-Free Condition 

Within a limited field of view we assume that the solar surface can be represented 
by the heliographic plane, which is tangent to the solar surface at the position 
given by the central meridian angle, the latitude, and the P- and Strangles. Over 
the heliographic plane, the image components of the magnetic field are related 
to the heliographic components according to the following transformation: 



B x — an-BJ, + ai 2 By + a\^B\ , 

By = a 2 i4+a22S;+ a 23Si> (1) 
B z = az\B l x + a 32 By + a 33 51, , 



where the coefficients are taken from Equation (1) of Gary and Hagyard 
(1990), B\ = jBii is the line-of-sight component of the field, and the components 
of the field perpendicular to the line-of-sight direction are given by 



B\ = Bx cos £ , and B y = B x sin £ , (2) 

where B^ is the magnitude of the component transverse to the line-of-sight and 
the azimuthal angle, £, is the angle between the transverse component of the field 
and the x'-axis (where x 1 and y 1 are coordinates on the plane perpendicular to 
the line-of-sight). The azimuthal angle can only be determined observationally 
within the range < £ < 180°. 

The derivative in the direction perpendicular to the heliographic plane, d/dz h , 
cannot be directly measured away from disk centre. However, the derivative 
along the line-of-sight direction, d/dz 1 , can be measured (see below for further 
discussion). These derivatives are related as follows: 



df df df df 

d? =ai3 d^ +a23 dy^ + a33 d^^ (3) 

where / is any differentiable function, z 1 is the distance along the line-of-sight 
direction, and x h and y h are coordinates on the heliographic plane. Derivatives 
with respect to horizontal heliographic distance, d/dx h and d/dy h , can be ap- 
proximated with techniques such as finite differences and finite elements using 
discrete measurements over the heliographic plane. For convenience, we work 
with the line-of-sight distance z 1 , although in practice the field is likely to be 
inferred as a function of optical depth; in which case a model atmosphere must be 
used to determine the correspondence between line-of-sight distance and optical 
depth (e.g., Maltby et ai, 1986; Vernazza, Avrett, and Loeser, 1981; Collados 
et ai, 1994; Socas-Navarro, 2007). 

Equations (1) and (3) can be used to express the divergence of the field in 
terms of observable quantities (i.e., derivatives of the image components of the 
field with respect to x h , y h , and z 1 ), 
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a 33 V B = D a + a 33 



dB\ 



(4) 



where 



D a = a 3 i—f- + a 32 -^-4- + (011033 - 013031) 



+ (021033 - 023031) — -f + (022033 - 023032) 





dy h ' 



+ (012033 - 013032) 



(5) 



Equations (4) and (5) show that the computation of the divergence requires 
the linc-of-sight derivatives of all three components of the magnetic field vector, 
except at disk centre where the coefficients a%j = Sij (see paper I for further 
discussion). 

Several techniques have appeared in the literature to measure the line-of-sight 
variation of the magnetic field (often but not always including the full magnetic 
field vector). We briefly discuss three that could be employed to evaluate the 
various terms in Equations (4) and (5). The first and possibly most straight- 
forward is to measure the field at different heights using multiple spectral lines 
(e.g., Liu et al, 1996; Leka and Metcalf, 2003). The inferred height difference 
between measurements can be several hundred kilometres or more depending on 
line choice, although exact heights cannot be easily quantified due to the lines' 
often broad range of formation height. A second approach makes use of a single 
(generally chromospheric) line that forms over a range of heights. Inversions 
performed at multiple, but limited wavelength ranges near (for example) the 
core and separately in the line wing can return the field strengths at different 
atmospheric heights (e.g., Metcalf et al., 1995; Eibe et al., 2002). Again, es- 
timates of the field are obtained at discrete but broad heights which may be 
separated by tens or hundreds of kilometres. A third approach recovers the line- 
of-sight stratification of the magnetic field during the inversion of the Stokes 
spectra through a formalism based on response functions (e.g., Ruiz Cobo and 
del Toro Iniesta, 1992; Collados et al., 1994; Del Toro Iniesta and Ruiz Cobo, 
1996; Wcstendorp Plaza et al., 1998, 2001; Socas-Navarro, Trujillo Bueno, and 
Ruiz Cobo, 2000; Socas-Navarro, 2005, 2007). The field is computed at discrete 
heights (inversion nodes) which may be separated by several hundred kilome- 
tres (e.g., Socas-Navarro, 2004); interpolation can map the results onto a grid 
with points separated by as little as 5 kilometres. For all of the aforementioned 
techniques, the field is effectively sampled discretely. Therefore, the line-of-sight 
derivatives in Equations (4) and (5) must be approximated. 

It is important to understand that the inferred line-of-sight derivatives of 
the transverse components of the field, dB^/dz 1 and dB^/dz 1 , depend on the 
ambiguity resolution at all heights used to approximate line-of-sight derivatives. 
Consequently, for ambiguity resolution algorithms based on these derivatives, 
vector magnetogram data should be disambiguated simultaneously at all relevant 
heights. 



cbl.tex; 3 November 2009; 23:43; p. 4 



Resolving the 180° Ambiguity in Vector Magnetograms 



5 



3. Synthetic Data and Performance Metrics 

Metcalf et al. (2006) compared several techniques for disambiguating single- 
height magnetogram data for which information about the variation of the field 
out of the heliographic plane was not available. To test the performance of the 
various algorithms based on the divergence we use the same magnetic field 
configurations that were used by Metcalf et al. (2006) for which the correct 
ambiguity resolution is known and magnetic field data are available at two 
heights. Methods based on the divergence require magnetic field data from at 
least two heights. All methods examined in this paper could be adapted for data 
with more than two heights; however, this can add a degree of complexity to the 
disambiguation process (depending on the algorithm) because all heights may 
need to be disambiguated simultaneously. For each synthetic magnetogram, the 
field is sampled at discrete spatial locations. The grid points are intended to 
mimic pixel centres, although the exact value of the field is provided. Thus, the 
effects of noise and unresolved structure in both the line-of-sight and horizontal 
heliographic directions are neglected (see e.g., Leka et al., 2009); these issues 
are beyond the scope of this investigation and will be considered in a separate 
paper. The grid points at each height are aligned in the line-of-sight direction 
and separated by a line-of-sight distance, Az 1 . 
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Figure 1. Synthetic vector magnetic field data. In both panels the underlying "continuum 
image" is a reverse-color image of B 2 ; positive/negative vertical magnetic flux is indicated 
by red/ blue contours at 100, 1000, 2000 G, and the magnetic neutral line is indicated by 
the black contour. Horizontal magnetic field is plotted at every fifth pixel, with magnitude 
proportional to arrow length. Tickmarks are in units of pixels, (a) Vector magnetic field at 
z = 0.00625L from the numerical simulation of Fan and Gibson (2004) at their time step 56 (L 
is a length scale, effectively the size of the simulation domain in the {/-direction). The pixel size 
is 0.00625L. This is an exact reproduction of Figure 1 from Metcalf et al. (2006). (b) Vector 
magnetic field for the multipole structure in image coordinates. The pixel size is 0.5". This is 
an exact reproduction of Figure 2 from Metcalf et al. (2006). 



The field shown in Figure 1(a) is a snapshot from the numerical MHD simu- 
lation of Fan and Gibson (2004) of a highly twisted flux tube emerging into an 
overlying potential field "arcade" . Because the simulation results arc provided 
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on a three-dimensional rectangular grid this case corresponds to disk centre, 
where the line-of-sight and z h -directions arc parallel. Figure 1(a) shows the field 
at the first height (as used by Metcalf et al., 2006); the second height used in 
this investigation lies above this, and is discussed in detail below. 

The field displayed in Figure 1(b) is constructed from a collection of point 
sources located on a plane below the surface. The contribution to the field from 
each source is calculated using the Green's function given by Chiu and Hilton 
(1977) with a constant value of the force-free parameter a. This magnctogram 
is located away from disk centre, centred at N 18°, W 45°, and it includes the 
effects of curvature. Again Figure 1(b) shows the field at the first height (as 
used by Metcalf et ai, 2006); the second height lies above this and is discussed 
further below. 

To quantify the performance of the various algorithms we will use the same 
four metrics presented in Metcalf et al. (2006): the fraction of pixels correctly 
disambiguated A4 area i the fraction of magnetic flux with the correct solution 
■Afflux, the fraction of strong horizontal field regions (£?h > 500 G) with the 
correct solution yV( h , and a measure of the discrepancy in the vertical current 
between the retrieved and correct solution, Mj z ; note that Mj z can have values 
outside the range [0, 1] unlike the other three metrics. Thus the results presented 
here can be directly compared to those obtained for the algorithms examined by 
Metcalf et al. (2006). 

4. Testing the Algorithms 

In this section we test several disambiguation methods all based on the divergence-^ 
free condition but which employ different assumptions and implementation de- 
tails. First, we examine two methods that have previously appeared in the 
literature: the Wu and Ai criterion (e.g., Wu and Ai, 1990; Cuperman, Li, 
and Semel, 1993; Li, Cuperman, and Semel, 1993; Li, Amari, and Fan, 2007) 
and the sequential minimisation method (e.g., Boulmczaoud and Amari, 1999). 
We find that both of these methods can produce significant errors for both 
synthetic data sets. Consequently, we present and test a third algorithm, the 
global minimisation method. 

4.1. The Wu and Ai Criterion 

One non-iterative approach that was originally proposed by Wu and Ai (1990) 
and tested on synthetic data by Cuperman, Li, and Semel (1993), Li, Cuper- 
man, and Semel (1993), and Li, Amari, and Fan (2007), involves expressing the 
divergence-free condition as an inequality. For any position on the solar disk 
(except the limb) the divergence-free condition can be written as 



Assuming that the various derivatives can be determined reliably, given the 
sign of dBl/dz 1 the ambiguity can in principle be resolved by choosing the 




(6) 
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azimuthal angle such that D a satisfies Equation (6). If dB\/dz l = or D a = 
then Equation (6) cannot be used to resolve the ambiguity. Consequently, for 
discrete observations Equation (6) cannot be used if \dB\jdz l \ < S(dB 1 z /dz 1 ) or if 
\D a \ < S(D a ), where S is an estimate of the measurement error in the respective 
derivatives. One advantage of using Equation (6) to resolve the ambiguity is 
that only the sign of dB l z /dz l is required, not its magnitude; note, however, that 
the sign and magnitude of dB^/dz 1 and dBy/dz 1 are required away from disk 
centre (see paper I). One limitation of using Equation (6) is that the inferred 
divergence of the field is not exactly zero (in general) when approximated from 
discrete observations, so the inequality is not exact (Cuperman, Li, and Semel, 
1993). 




50 100 150 200 50 100 150 200 



Figure 2 . (a) Disambiguation results for the Wu and Ai criterion applied to the flux tube and 
arcade configuration (Figure 1(a)). Contours are the same as Figure 1(a), except that the green 
contour is the magnetic neutral line. Areas with the correct /incorrect ambiguity resolution are 
black/white (the same convention as Metcalf et ai, 2006). Initial configuration has B y > 0. 
(b) Same as (a) except the initial configuration of azimuthal angles is random. 

The algorithm proceeds as follows. At a given height each column of pixels 
is scanned from j = 1 to j = n y , and different columns are disambiguated from 
i — 1 to i — n x . Hereafter, the indices i, j, and k correspond to pixel labelling 
in the x 1 -, y 1 -, and z'-dircctions, respectively, and n x and n y refer to the number 
of pixels in the x 1 - and y'-directions, respectively. For 1 < i < n x and 1 < j < 
n y at a given height k we approximate the horizontal hcliographic derivatives 
d/dx h and d/dy h at pixel (i,j,k) with three-point finite differences using the 
measurements at pixels: k), k), and + k). At disk centre, where 

the heliographic and image coordinates are equivalent, these approximations 
reduce to standard two-point forward finite differences. At the boundaries i = n x 
and j = n y the geometry of the finite differencing stencil is modified to use only 
pixels within the field of view. At the lower/upper height forward/backward 
finite differences arc used to approximate the line-of-sight derivative d/dz 1 using 
measurements from the two heights. 

Away from disk centre the lower height k = 1 is examined first followed by 
the second height k = 2. At disk centre only the lower height is disambiguated; 
this is possible because dB x jdz i and dBy/dz 1 are not required to compute the 
divergence at disk centre (see Equations (4) and (5)) and the approximations for 
horizontal heliographic derivatives involve measurements at a single height. At 
each pixel (i, j, k) we evaluate the quantities in Equation (6). If the criterion is 
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Table 1. Performance metrics for the flux tube and arcade at disk centre for 
the various ambiguity resolution algorithms (with Az 1 = 1 pixel). The lower 
height is evaluated. 



Method 


A'farca 


A4 flux 


M h 


Mj z 


Wu and Ai, initialised with B l y > 


0.98 


0.98 


0.97 


-0.29 


Wu and Ai, initialised by randomising {; 


0.62 


0.63 


0.63 


-15.24 


Sequential minimisation 


0.94 


0.93 


0.93 


-0.88 


Global minimisation 


1.00 


1.00 


1.00 


1.00 



satisfied then no action is taken and we move to the next pixel in the sequence. 
If the criterion is violated then the other choice of azimuthal angle is selected 
at pixel (i, j, k). To be consistent with the algorithms described in Wu and Ai 
(1990), Cuperman, Li, and Scmel (1993), Li, Cuperman, and Semel (1993), and 
Li, Amari, and Fan (2007) we do not check that IdB^/dz 1 ] and \D a \ exceed their 
respective measurement errors nor that Equation (6) is satisfied if the other 
choice of azimuthal angle is selected. 

To demonstrate the sensitivity of this algorithm to the initial configuration of 
azimuthal angles we show two cases. First, we initialise the data by setting B y > 
at each pixel; this effectively smoothes the initial configuration. At disk centre 
this is essentially the approach taken by Li, Amari, and Fan (2007), who also used 
synthetic data from the numerical simulation of Fan and Gibson (2004), although 
there may be some subtle differences between our implementation and the Li, 
Amari, and Fan implementation (such as the order in which pixels are examined 
and the exact heights used). Second, we initialise the data by randomising the 
choice of azimuthal angle at each pixel (the initial choice of azimuthal angle is 
changed in approximately 50% of pixels) . The performance of this algorithm for 
both initialisations is shown in Figures 2 and 3, and Tables 1 and 2. 

Evidently, this method is very sensitive to the smoothness of the initial con- 
figuration. This is because the approximation for D a may be misleading when 
the initial configuration is not smooth. This algorithm is also sensitive to the 
order in which pixels are examined: for the off-disk-ccntre data (Figure 3) the 
results at the upper height tend to be worse than those at the lower height. This 
is because the smoothness of the initial configuration can be perturbed during 
the disambiguation process, resulting in a misleading approximation for D a . 

The assumption made by this method is tested by using the correct solution 
to evaluate each of the terms in Equation (6) with the various derivatives approx- 
imated as described above. In Figures 4 and 5 areas that satisfy the inequality 
in Equation (6) for the correct solution are shown along with areas that violate 
Equation (6). Also shown are areas which suffer from the degenerate case where 
Equation (6) cannot be used to distinguish the two realisations, as the sign of 
D a is the same for both the correct and the incorrect choice of azimuthal angle 
at the pixel of interest. This is the case for a large fraction of the field of view 
for both synthetic data sets. The method appears to do quite well in Figure 2(a) 
and Figure 3(a) because we do not test the other realisation if the initial one 
satisfies Equation (6) and we do not check that Equation (6) is satisfied after the 
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Figure 3. (a) Disambiguation results for the Wu and Ai criterion applied to the multipolc 
field configuration (Figure 1(b)), following Figure 2(a). (b) Same as (a) except for the second 
height, 1 pixel (or 0.5") above that shown in (a), (c) Same as (a) except the initial configuration 
of azimuthal angles is random, (d) Same as (c) except for the second height. 

choice of azimuthal angle is changed. The degeneracy is partly due to the details 
of the implementation and could be avoided by changing the azimuthal angle 
in all pixels referenced by the finite differencing stencil at a particular location, 
rather than at a single pixel; unfortunately, such an approach would be very 
sensitive to the initial configuration. 

4.2. The Sequential Minimisation Method 

Boulmezaoud and Amari (1999) presented a method which assumes that the 
magnitude of the divergence is minimised over the set of pixels used to approx- 
imate the divergence at a particular location. Disambiguation then proceeds as 
a sequence of minimisation problems based at each pixel in the field of view. 

We apply this approach as follows. At a given height, we examine pixels in the 
same order as in Section 4.1. For the pixels referenced by the finite differencing 
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Table 2. Performance metrics for the multipole field positioned away from 
disk centre for the various ambiguity resolution algorithms (Az 1 = 1 pixel). The 
lower height is evaluated. 



Method 


A'farca 


A4 flux 


M h 


Mj, 


Wu and Ai, initialised with B y > 


0.92 


0.94 


0.93 


-1.14 


Wu and Ai, initialised by randomising £ 


0.61 


0.59 


0.66 


-13.70 


Sequential minimisation 


0.95 


0.93 


0.94 


-3.36 


Global minimisation 


1.00 


1.00 


1.00 


1.00 




Figure 4. The Wu and Ai criterion applied to the correct solution for the flux tube and arcade 
configuration (Figure 1(a)). Contours are the same as Figure 2. Areas that satisfy/ violate the 
Wu and Ai criterion, Equation (6), are coloured black/ white. Areas for which the sign of D a 
is the same for both choices of the azimuthal angle are coloured gray. 

stencil at a particular location we choose the permutation of azimuthal angles 
that corresponds to the minimum of KV-B)^-^!, the approximation for the 
divergence at the pixel (i,j,k). We compute the divergence using the methods 
described in Section 4.1, with all approximations made about the lower height. 
The azimuthal angle at each pixel is fixed by the stencil that first references it. 
This approach examines both permutations of the azimuthal angle at each pixel, 
therefore the initial configuration does not affect its performance, although it is 
sensitive to the order in which pixels are examined. 

The performance of this method is shown for the flux tube and arcade (Fig- 
ure 6(a)). At disk centre the disambiguation at each height is independent 
because dB^/dz 1 and dB y jdz l arc not required to compute the divergence and 
the approximations for horizontal hcliographic derivatives involve measurements 
at a single height; results for the lower height arc displayed in Figure 6(a). For 
the multipole field positioned away from disk centre (Figures 7(a) and (b)) both 
heights are disambiguated simultaneously as required for off-disk-centre data. 
Performance metrics for the lower height in each case are provided in Tables 1 
and 2. 

For each pixel the ambiguity is resolved by the stencil that references it first. 
Hence, for the pixels in a given stencil this algorithm assumes that the correct 
ambiguity resolution will be recovered for the unresolved pixels if the previously 
resolved pixels have the correct ambiguity resolution. We show areas for which 



cbl.tex; 3 November 2009; 23:43; p. 10 



Resolving the 180° Ambiguity in Vector Magnetograms 



11 




50 100 150 200 250 300 50 100 150 200 250 300 



Figure 5. (a) The Wu and Ai criterion applied to the correct solution for multipole field 
configuration at the lower height (Figure 1(b)); black/ white/ gray indicate areas which satisfy, 
do not satisfy, and are degenerate to the inequality in Equation (6), following Figure 4. (b) 
Same as (a) except for the second height. 




Figure 6. (a) Disambiguation results for the sequential minimisation method applied to the 
flux tube and arcade configuration, following Figure 2. (b) Validity of the assumption made by 
the sequential minimisation method. Areas that satisfy /violate the assumption are coloured 
black/ white. Contours are the same as (a). 



this assumption is valid for both test cases in Figure 6(b) and Figures 7(c) and 
(d). Evidently, this assumption is not universally valid when the divergence is 
approximated from discrete observations. However, the fraction of pixels that 
violate this assumption is smaller than the fraction that retrieve the incorrect 
azimuthal angle. There arc two reasons that this method fails at a particular 
pixel: i) the assumption is incorrect; and ii) the propagation of an erroneous 
solution that originates in a region where the assumption is incorrect. Prop- 
agation accounts for most of the areas with the incorrect azimuthal angle in 
Figure 6(a) and Figures 7(a) and (b). Most, but not all, pixels that violate 
the assumption initiate the propagation of errors. Also, there are pixels that 
violate the assumption but still retrieve the correct ambiguity resolution; these 
situations can occur when the assumption is violated and a neighbouring pixel 
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Figure 7. (a) Disambiguation results for the sequential minimisation method applied to the 
multipolc field configuration, following Figure 3(a). (b) Same as (a) except for the second 
height, 1 pixel (or 0.5") above that shown in (a), (c) Validity of the assumption made by the 
sequential minimisation method applied to the multipolc field, following Figure 6(b). (d) Same 
as (c) except for the second height. 

has the incorrect ambiguity resolution, allowing the algorithm to get back onto 
the correct solution. 

The algorithm described above is different from the implementation of Boul- 
mezaoud and Amari (1999) in that they used four-point finite elements to 
approximate horizontal hcliographic derivatives (we note that there may also 
be other differences) . We have tested other versions of this algorithm employing 
different approximations for the horizontal hcliographic derivatives, such as four- 
point finite differences and finite elements. The results are very similar to those 
shown in Figure 6 and 7: all versions examined violate the underlying assumption 
at some locations and subsequently propagate errors in the solution. 

This method can propagate errors because it is sequential. We have tried a 
non-sequential approach which assumes that (V-B)^-^! is minimised within the 
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stencil based at each pixel regardless of the solution obtained in neighbouring 
pixels. For both synthetic data sets the fraction of pixels for which this assump- 
tion is violated is actually larger than for the sequential approach. Therefore, we 
conclude that the non-sequential approach does not substantially improve the 
performance of this method. 

4.3. The Global Minimisation Method 

In this section we present an iterative, global method which is based on the 
approximation for the divergence summed over the entire field of view. This 
approach is related to the "minimum energy" algorithm for single-height data 
(see e.g., Metcalf, 1994; Metcalf et al., 2006; Leka et al., 2009; Leka, Barnes, and 
Crouch, 2009). We assume that the correct configuration of azimuthal angles 
over the field of view is the one that corresponds to the minimum of 

^ = EE(i( v - B )^i + i( v - B )^+ii) ■ (7) 

»=i j=i 

Equation (7) consists of two terms corresponding to the approximations for the 
divergence based at each of the two heights, k and k + 1. For data positioned 
away from disk centre, where both heights are disambiguated simultaneously, we 
find that algorithms based on Equation (7) consistently retrieve a better solution 
than algorithms based on a single term. For the divergence at the first/second 
height the line-of-sight derivatives are approximated with forward/backward 
finite differences. For each of these terms the horizontal heliographic derivatives 
are approximated using three-point finite differences as described in Section 4.1. 
We have tested several different ways to approximate horizontal heliographic 
derivatives, such as four-point finite differences and three- and four-point finite 
elements, but these do not significantly affect the performance of the algorithm. 

The number of possible solutions for this optimisation problem is finite but 
very large, 2 n ^ nv , making a brute-force search for the global minimum of E 
impractical. One method that works well for large, discrete, global optimisation 
problems such as this is simulated annealing {e.g., Metropolis et al., 1953; 
Kirkpatrick et al., 1983; Press et al., 1992). Using simulated annealing, the 
global minimisation method proceeds as follows. A pixel is selected at random. 
If changing the azimuthal angle at that pixel reduces the value of E then the 
change is always accepted. If changing the azimuthal angle increases the value 
of E then the change is sometimes accepted, with a probability given by 

p = cxp (-AE/T) , (8) 

where T is the "temperature" and AE is the change in E (Equation (7)) 
caused by changing the azimuthal angle; positive values of AE imply that E 
increases due to the reconfiguration. The ability to sometimes accept an "uphill" 
reconfiguration enables simulated annealing to avoid getting stuck in a sub- 
optimal solution, provided that the cooling rate is sufficiently slow and the initial 
temperature is sufficiently high (e.g., Geman and Geman, 1984). 
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The annealing schedule (the initial temperature, the number of reconfigura- 
tions attempted at each temperature setting, and the cooling rate) can require 
some experimentation; in the interests of simplicity we use the following pre- 
scription. The initial temperature, To, should greatly exceed any expected value 
of AE. Starting with the initial configuration of azimuthal angles we examine 
I00n x n y reconfigurations of randomly chosen pixels (which are all accepted) 
and take two times the maximum value of \AE\ as the initial temperature. It 
should be noted that this process effectively randomises the initial configuration 
of the azimuthal angles and therefore the initial configuration does not affect the 
final result (likewise, the high initial temperature further randomises the initial 
configuration because almost all reconfigurations will be accepted during early 
temperature increments). At each temperature setting we examine 20n x n y recon- 
figurations of the azimuthal angle at randomly chosen pixels. The temperature 
varies according to T t = c l To, where c is the cooling rate and t is the label of the 
temperature increment; we use c = 0.999. We continue in this fashion until one 
of the following three conditions is satisfied: i) no reconfigurations are accepted 
at the most recent temperature setting; ii) the temperature is sufficiently low 
(less than 10 -7 To); or Hi) for ten consecutive temperature increments the value 
of E does not change substantially (i.e., \E t — E t -i\/ (\E t \ + l-Et-ii) < 10~ 5 , 
where E t is the value of E at the end of temperature increment t). 

For the flux tube and arcade positioned at disk centre, we show the fraction 
of pixels correctly disambiguated by the global minimisation method, A^areaj as 
a function of the vertical grid spacing Az 1 ( Figure 8); only the lower height 
is evaluated because the disambiguations at each height are independent at 
disk centre for the derivative approximations employed here. The simulation 
results of Fan and Gibson (2004) have a vertical grid spacing Az 1 of 0.00625L 
(equal to the horizontal grid spacing). For the first height we use the snapshot 
displayed in Figure 1(a) at z 1 = 0.00625L. We allow the second height to vary, 
using all heights in the range z 1 = 0.0125L, . . . , 0.625L. For each setting of Az 1 
we disambiguate twenty times with the algorithm described above but with a 
different sequence of random numbers for each disambiguation. 

Independent disambiguations can behave differently because simulated an- 
nealing is a stochastic optimisation algorithm: the decision to move to higher 
E is probabilistic (Equation (8)) and the order in which pixels are examined is 
random. Consequently, if the annealing schedule is not ideal (e.g., if the cooling 
rate is too fast) the global minimum may not be found and a different answer may 
be obtained for a different sequence of random numbers. To ensure that the global 
minimum is found for real magnetogram data it may be necessary to repeat 
the disambiguation process multiple times with different sequences of random 
numbers, a slower cooling rate, and/or with more attempted reconfigurations at 
each temperature setting. 

For the flux tube and arcade at disk centre with Az 1 = 1 pixel the global 
minimisation method retrieves the correct azimuthal angle at every pixel for 
all twenty disambiguations with the imposed annealing schedule (see Table 1). 
Several points are worth noting regarding Figure 8: i) Broadly speaking, the 
performance decreases as Az 1 increases, ii) The scatter of the A^arca values 
increases as Az 1 increases (most noticeably for Az 1 > 60 pixels). Hi) For Az 1 > 
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Figure 8. The fraction of pixels correctly disambiguated at the lower height, A4 area j as a 
function of vertical grid spacing Az 1 (in units of horizontal pixels) for the global minimisation 
method with a fixed annealing schedule applied to the flux tube and arcade configuration 
(Figure 1(a)). For each setting of Az 1 twenty disambiguations are executed with different 
sequences of random numbers. The gray dots represent the resultant .Marea value for each of 
these disambiguations. The black line connects the median of the .Marea values for each Az 1 
setting and the gray lines connect the corresponding minimum and maximum of the .Marea 
values. The dotted horizontal line at .Marea = 1 is included for reference. 



3 pixels the correct solution is no longer found for all twenty disambiguations. 
The maximum of the yV/f aroa values is one at this point so the correct solution 
could be obtained simply by disambiguating multiple times and selecting the 
solution with the smallest value of E. iv) The point where the maximum of the 
.Marea values deviates from one is Az 1 10 pixels. Beyond this point, a perfectly 
correct solution is not obtained for any of the twenty disambiguations with the 
imposed annealing schedule, yet for the largest values of Az 1 the median value 
for .Marea is still greater than 0.99. 




50 100 150 200 50 100 150 200 

Figure 9. Disambiguation results for the global minimisation method applied to the flux tube 
and arcade configuration with Az 1 = 65 pixels, following Figure 2. (a) The solution retrieved 
with the best .Marea value, (b) The solution retrieved with the worst .Marea value. 



For large values of Az 1 the difference between the best and worst solutions 
retrieved by the global minimisation method can be significant (Figure 9). The 
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main cause of such deviations is the annealing schedule. We have conducted 
experiments that show that more consistent and better results can be obtained 
at large values of Az 1 with annealing schedules that involve slower cooling rates 
or more attempted reconfigurations at each temperature setting. 

For the multipole field positioned away from disk centre the variation of Ai area 
as a function of the line-of-sight grid spacing Az 1 is shown (Figure 10); both 
heights are disambiguated simultaneously as required for off-disk-centre data. 
For the first height we use the field displayed in Figure 1(b). We allow the second 
height to vary in increments along the line-of-sight direction corresponding to 
the horizontal pixel size, 0.5". Of the twenty disambiguations with Az 1 = 1 pixel, 
only two solutions obtained A4 a rca = 1 for both heights (see Table 2); the 
remainder of solutions obtained .Marca = 0.99 or better at each height. 

The performance of the global minimisation method applied to the multipole 
field configuration (Figure 10) degrades more rapidly as Az 1 increases than it 
does for the flux tube and arcade (Figure 8). This is partly because the op- 
timisation problem for data positioned away from disk centre is substantially 
more challenging, as both heights must be disambiguated simultaneously. Dif- 
ferences in the variation of the two fields with height may also contribute to the 
contrasting performance. 

The synthetic data for the multipole field configuration include the effects of 
curvature but we have not included curvature in the treatment of the divergence 
(Equations (4) and (5)). Experiments with two synthetic data sets that have the 
same field configuration but include and exclude curvature indicate that curva- 
ture does not significantly influence the performance of the global minimisation 
algorithm, provided that the field of view is much smaller than the radius of 
curvature of the Sun. 

In Figure 10, we see that for some values of Az 1 larger than 31 pixels (red 
dots in Figure 10), the algorithm retrieves solutions with E less than that for 
the correct solution. Thus, the assumption that the correct solution corresponds 
to the configuration of azimuthal angles that minimises E (Equation (7)) is 
violated in these cases. Nevertheless, the performance of the resulting solutions 
is quite good, with M ave a, values greater than 0.99 at both heights (see e.g., 
Figure 11). For Az 1 larger than 53 pixels no such solutions arc found, which 
probably indicates that the fixed annealing schedule employed here is less than 
ideal for large values of Az 1 . 

Results for a case where the retrieved solution has E less than that for the 
correct solution are shown for the multipole configuration with Az 1 = 50 pixels 
(Figure 11). Interestingly areas that fail occur at different locations at each 
height; this situation is fairly typical for Az 1 > 25 pixels. For smaller values of 
Az 1 incorrect regions tend to occur at similar locations at both heights. 

For the flux tube and arcade we have not found a solution for which E is less 
than that for the correct solution, in the case where the horizontal resolution is 
fixed but the vertical resolution varies (Figure 8). We find that this situation does 
arise for the flux tube and arcade when both the horizontal and vertical sampling 
are about a factor of seven times coarser than on the original grid. Clearly, the 
validity of the underlying assumption made by the global minimisation method 
- that the correct configuration of azimuthal angles corresponds to the minimum 
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Fi gure 10. (a) The fraction of pixels correctly disambiguated at the lower height, .A/Iarca, 
as a function of vertical grid spacing A2 1 (in units of horizontal pixels) for the global min- 
imisation method with a fixed annealing schedule applied to the multipole field configuration 
in Figure 1(b), following Figure 8. The red dots represent the jVfarea values for each of the 
disambiguations that retrieved a solution with E less than that for the correct solution, (b) 
Same as (a) except for the second height. 
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Figure 11. (a) Disambiguation results for the global minimisation method applied to the 
multipole field configuration with A2 1 = 50 pixels, following Figure 3(a). For this particular 
case E for the retrieved solution is less than that for the correct solution, (b) Same as (a) 
except for the second height, 50 pixels above that shown in (a). 



value of E - is sensitive to the spatial resolution of the observations. This sug- 
gests that for real magnetogram data (that may contain noise and/or unresolved 
structure) additional constraints or criteria may be required to determine if the 
global minimum corresponds to the correct solution. It may also be possible 
to use measurement error information to determine the uncertainties associated 
with the configuration corresponding to the global minimum. 



cbl.tex; 3 November 2009; 23:43; p. 17 



18 



A.D. Crouch et al. 



5. Conclusions 

We investigate how the divergence-free property of magnetic fields can be ex- 
ploited to resolve the azimuthal ambiguity that is present in solar vector mag- 
nctogram data by using linc-of-sight and horizontal hcliographic derivative in- 
formation approximated from discrete measurements. Using synthetic data we 
test several methods that each make different assumptions about how to use 
the divergence to resolve the ambiguity. We find that all of the assumptions 
considered can be violated when the divergence is approximated from discrete 
measurements (depending on the nature of the field and the resolution of the 
observations), but some assumptions are more robust than others. 

The Wu and Ai criterion expresses the divergence-free condition as an in- 
equality {e.g., Wu and Ai, 1990; Cuperman, Li, and Semel, 1993; Li, Cuperman, 
and Semel, 1993; Li, Amari, and Fan, 2007). We find that methods based on 
this criterion are very sensitive to the smoothness of the initial configuration of 
azimuthal angles and the order in which pixels are examined. In addition, the 
underlying assumption can be incorrect over a significant fraction of the field 
of view when derivatives are approximated from discrete measurements. The 
sequential minimisation method (Boulmczaoud and Amari, 1999) assumes that 
the magnitude of the divergence is minimised over the pixels used to approximate 
the divergence at each point. Our tests indicate that this approach is more robust 
than the Wu and Ai criterion but can produce errors when the underlying as- 
sumption is violated. Moreover, the algorithm can propagate erroneous solutions 
into regions that do not violate the underlying assumption. 

We find that the most promising algorithm is the global minimisation method 
which assumes that the correct configuration of azimuthal angles corresponds to 
the minimum of the magnitude of the divergence summed over the entire field 
of view. This method relies on the magnitude of the divergence. Therefore, with 
this method, the sign and magnitude of the line-of-sight derivatives of all three 
components of the magnetic field vector are required to resolve the azimuthal 
ambiguity away from disk centre. 

For the synthetic data employed here, the field is sampled at discrete spatial 
locations where the exact value of the field is provided. Thus, the effects of 
noise and unresolved structure in both the horizontal hcliographic and line- 
of-sight directions are neglected. The performance of the global minimisation 
method is very promising, but before this approach can be reliably applied to 
real magnetogram data its performance should be tested with more realistic 
synthetic data that include the effects of noise and unresolved structure {e.g., 
Leka et al, 2009). 
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